##Copy the part starting with source() into the console and run it there - this will create a log file that includes the output

#source("~/work/November_2024/Script_MainPlots.R", echo = T, max.deparse.length=10000)
rm(list=ls(all=TRUE))
sink("~/work/November_2024/R_out/log_files/full_log_mainPLOTS25.log", split = T)

# ======================================================================
# Project:    Closed borders, closed minds? COVID-related border closures,
#             EU support and hostility towards immigrants
#
# Script:     Figure 3 + Figure 4
#
# Authors:    Lisa Herbig (l.j.herbig@uva.nl)
#             Asli Unan (a.unan@uva.nl)
#
# Date:       2025-03-24
# ======================================================================

# Description:
# This script generates Figure 3 and Figure 4, illustrating the average treatment effect on attitudes towards refugees and attachment with Europe. 

# ----------------------------------------------------------------------
# Load libraries
# ----------------------------------------------------------------------
library("did2s")
library("dplyr")
library("foreign")
library("ggplot2")
library("haven")
library("lmtest")
library("lubridate")
library("modelsummary")
library("PanelMatch")
library("panelView")
library("readr")    
library("tidyr")
library("tidyverse")

#############################################START SCRIPT###############################################################
## Read the files in
df_iso_work_dd = readRDS("~/work/November_2024/R_out/df_iso_work_dd_25.rds")


## Define output folder
output_folder <- "~/work/November_2024/R_out/main_plots/"


#############################################CONNECTION EUROPE
# No Matching
any_closure_binary_control0_nm = PanelMatch(lag = 4, time.id = "weekyear_i",
                                            unit.id = "iso2.y_i", 
                                            treatment = "any_closure_binary_control0", 
                                            refinement.method = "none", 
                                            data = as.data.frame(df_iso_work_dd), qoi = "att" ,
                                            covs.formula = ~ East + age_n + female_n + unemployed_n + working_abroad_n + university_n + EU_residence_n  + kr_foreigner_n + CONTACTABROAD_n + HOMEOFFICE_n + Voting_2018_LeftRight_n+ mode_phone_n +  iso2_pop_n + I(lag(kr_inz_rate_n, 1:4)),
                                            outcome.var = "con_europe_n",
                                            size.match = 1,
                                            match.missing = FALSE,
                                            matching = FALSE,
                                            lead = 1:4, forbid.treatment.reversal = FALSE, 
                                            placebo.test = T)

pt_DL_NN <- placebo_test(any_closure_binary_control0_nm, data = as.data.frame(df_iso_work_dd), plot = F, se.method = "conditional")
pt_DL_NN_t_4_est <- pt_DL_NN$estimates[1]
pt_DL_NN_t_3_est <- pt_DL_NN$estimates[2]
pt_DL_NN_t_2_est <- pt_DL_NN$estimates[3]

pt_DL_NN_t_4_high <- pt_DL_NN_t_4_est+(pt_DL_NN$standard.errors[1]*1.96)
pt_DL_NN_t_3_high <- pt_DL_NN_t_3_est+(pt_DL_NN$standard.errors[2]*1.96)
pt_DL_NN_t_2_high <- pt_DL_NN_t_2_est+(pt_DL_NN$standard.errors[3]*1.96)
pt_DL_NN_t_4_low <- pt_DL_NN_t_4_est-(pt_DL_NN$standard.errors[1]*1.96)
pt_DL_NN_t_3_low <- pt_DL_NN_t_3_est-(pt_DL_NN$standard.errors[2]*1.96)
pt_DL_NN_t_2_low <- pt_DL_NN_t_2_est-(pt_DL_NN$standard.errors[3]*1.96)

Results_DL_NN <- PanelEstimate(sets = any_closure_binary_control0_nm, data = as.data.frame(df_iso_work_dd), se.method = "conditional")
summary(Results_DL_NN)

# Matching
any_closure_binary_control0_m = PanelMatch(lag = 4, time.id = "weekyear_i",
                                           unit.id = "iso2.y_i", 
                                           treatment = "any_closure_binary_control0", 
                                           refinement.method = "mahalanobis", 
                                           data = as.data.frame(df_iso_work_dd), qoi = "att",
                                           covs.formula = ~ East + age_n + female_n + unemployed_n + working_abroad_n + university_n + EU_residence_n  + kr_foreigner_n + CONTACTABROAD_n + HOMEOFFICE_n + Voting_2018_LeftRight_n + mode_phone_n +  iso2_pop_n + I(lag(kr_inz_rate_n, 1:4)),
                                           outcome.var = "con_europe_n",
                                           size.match = 1,
                                           match.missing = TRUE,
                                           matching = TRUE,
                                           lead = 1:4, forbid.treatment.reversal = FALSE, 
                                           placebo.test = T)


pt_DL_mahalanobis <- placebo_test(any_closure_binary_control0_m, data = as.data.frame(df_iso_work_dd), plot = F)

pt_DL_mahalanobis_t_4_est <- pt_DL_mahalanobis$estimates[1]
pt_DL_mahalanobis_t_3_est <- pt_DL_mahalanobis$estimates[2]
pt_DL_mahalanobis_t_2_est <- pt_DL_mahalanobis$estimates[3]

pt_DL_mahalanobis_t_4_high <- pt_DL_mahalanobis_t_4_est+(pt_DL_mahalanobis$standard.errors[1]*1.96)
pt_DL_mahalanobis_t_3_high <- pt_DL_mahalanobis_t_3_est+(pt_DL_mahalanobis$standard.errors[2]*1.96)
pt_DL_mahalanobis_t_2_high <- pt_DL_mahalanobis_t_2_est+(pt_DL_mahalanobis$standard.errors[3]*1.96)
pt_DL_mahalanobis_t_4_low <- pt_DL_mahalanobis_t_4_est-(pt_DL_mahalanobis$standard.errors[1]*1.96)
pt_DL_mahalanobis_t_3_low <- pt_DL_mahalanobis_t_3_est-(pt_DL_mahalanobis$standard.errors[2]*1.96)
pt_DL_mahalanobis_t_2_low <- pt_DL_mahalanobis_t_2_est-(pt_DL_mahalanobis$standard.errors[3]*1.96)

Results_DL_mahalanobis <- PanelEstimate(sets = any_closure_binary_control0_m, data = as.data.frame(df_iso_work_dd), se.method = "conditional")
summary(Results_DL_mahalanobis)


# Generate Figure 3

plot_table_NN <- rbind(pt_DL_NN_t_4_est, pt_DL_NN_t_3_est, pt_DL_NN_t_2_est, 0, as.data.frame(Results_DL_NN$estimates))
colnames(plot_table_NN) <- "estimate"
plot_table_NN$high <- c(pt_DL_NN_t_4_high, pt_DL_NN_t_3_high, pt_DL_NN_t_2_high, 0, Results_DL_NN$estimates + (Results_DL_NN$standard.error*1.96))
plot_table_NN$low <- c(pt_DL_NN_t_4_low, pt_DL_NN_t_3_low, pt_DL_NN_t_2_low, 0, Results_DL_NN$estimates - (Results_DL_NN$standard.error*1.96))
plot_table_NN$Matching <- "No Matching" 
rownames(plot_table_NN)[1:5] <- c("t-4","t-3","t-2","t-1","t")
plot_table_NN$Week <- rownames(plot_table_NN)


plot_table_R <- rbind(pt_DL_mahalanobis_t_4_est, pt_DL_mahalanobis_t_3_est, pt_DL_mahalanobis_t_2_est, 0, as.data.frame(Results_DL_mahalanobis$estimates))
colnames(plot_table_R) <- "estimate"
plot_table_R$high <- c(pt_DL_mahalanobis_t_4_high, pt_DL_mahalanobis_t_3_high, pt_DL_mahalanobis_t_2_high, 0, Results_DL_mahalanobis$estimates + (Results_DL_mahalanobis$standard.error*1.96))
plot_table_R$low <- c(pt_DL_mahalanobis_t_4_low, pt_DL_mahalanobis_t_3_low, pt_DL_mahalanobis_t_2_low, 0, Results_DL_mahalanobis$estimates - (Results_DL_mahalanobis$standard.error*1.96))
plot_table_R$Matching <- "Matching + Refinement" 
rownames(plot_table_R)[1:5] <- c("t-4","t-3","t-2","t-1","t")
plot_table_R$Week <- rownames(plot_table_R)

plot_table_DL <- rbind(plot_table_R, plot_table_NN)
plot_table_DL$Week <- factor(plot_table_DL$Week, levels = c("t-4", "t-3", "t-2", "t-1", "t", "t+1", "t+2", "t+3"))

plot_table_DL$`Estimation Strategy` <- factor(plot_table_DL$Matching, levels = c("No Matching", "Matching + Refinement"))


# Plot
pd = position_dodge(.6)   
ggplot(plot_table_DL, aes(x = Week, y = estimate, color = `Estimation Strategy`)) + geom_hline(yintercept = 0, color = "gray40", linetype="dashed") +
  geom_point(shape = 16, size  = 6, position = pd) +
  geom_errorbar(aes(ymin  = low, ymax  = high), width = 0, size  = 1.45, position = pd) +
  theme_bw() + theme() + ylab("Estimated Treatment Effect on Attachment with Europe") + scale_color_manual(values= c("gray70", "#377EB8")) + theme(legend.text=element_text(size=20), axis.text = element_text(size = 20), axis.title = element_text(size = 20), text = element_text(size = 20), plot.title = element_text(hjust = 0.5)) + ggtitle("Border Closures") 

ggsave(path = "~/work/November_2024/R_out/main_plots/", filename = "Con-Europe_MainPlot.pdf")




#############################################Refugee
# No Matching
any_closure_binary_control0_nm = PanelMatch(lag = 4, time.id = "weekyear_i",
                                            unit.id = "iso2.y_i", 
                                            treatment = "any_closure_binary_control0", 
                                            refinement.method = "none", 
                                            data = as.data.frame(df_iso_work_dd), qoi = "att" ,
                                            covs.formula = ~ East + age_n + female_n + unemployed_n + working_abroad_n + university_n + EU_residence_n  + kr_foreigner_n + CONTACTABROAD_n + HOMEOFFICE_n + Voting_2018_LeftRight_n+ mode_phone_n +  iso2_pop_n + I(lag(kr_inz_rate_n, 1:4)),
                                            outcome.var = "immigration_pca_n",
                                            size.match = 1,
                                            match.missing = FALSE,
                                            matching = FALSE,
                                            lead = 1:4, forbid.treatment.reversal = FALSE, 
                                            placebo.test = T)

pt_DL_NN <- placebo_test(any_closure_binary_control0_nm, data = as.data.frame(df_iso_work_dd), plot = F, se.method = "conditional")
pt_DL_NN_t_4_est <- pt_DL_NN$estimates[1]
pt_DL_NN_t_3_est <- pt_DL_NN$estimates[2]
pt_DL_NN_t_2_est <- pt_DL_NN$estimates[3]

pt_DL_NN_t_4_high <- pt_DL_NN_t_4_est+(pt_DL_NN$standard.errors[1]*1.96)
pt_DL_NN_t_3_high <- pt_DL_NN_t_3_est+(pt_DL_NN$standard.errors[2]*1.96)
pt_DL_NN_t_2_high <- pt_DL_NN_t_2_est+(pt_DL_NN$standard.errors[3]*1.96)
pt_DL_NN_t_4_low <- pt_DL_NN_t_4_est-(pt_DL_NN$standard.errors[1]*1.96)
pt_DL_NN_t_3_low <- pt_DL_NN_t_3_est-(pt_DL_NN$standard.errors[2]*1.96)
pt_DL_NN_t_2_low <- pt_DL_NN_t_2_est-(pt_DL_NN$standard.errors[3]*1.96)

Results_DL_NN <- PanelEstimate(sets = any_closure_binary_control0_nm, data = as.data.frame(df_iso_work_dd), se.method = "conditional")
summary(Results_DL_NN)

# Matching
any_closure_binary_control0_m = PanelMatch(lag = 4, time.id = "weekyear_i",
                                           unit.id = "iso2.y_i", 
                                           treatment = "any_closure_binary_control0", 
                                           refinement.method = "mahalanobis", 
                                           data = as.data.frame(df_iso_work_dd), qoi = "att",
                                           covs.formula = ~ East + age_n + female_n + unemployed_n + working_abroad_n + university_n + EU_residence_n  + kr_foreigner_n + CONTACTABROAD_n + HOMEOFFICE_n + Voting_2018_LeftRight_n + mode_phone_n +  iso2_pop_n + I(lag(kr_inz_rate_n, 1:4)),
                                           outcome.var = "immigration_pca_n",
                                           size.match = 1,
                                           match.missing = TRUE,
                                           matching = TRUE,
                                           lead = 1:4, forbid.treatment.reversal = FALSE, 
                                           placebo.test = T)


pt_DL_mahalanobis <- placebo_test(any_closure_binary_control0_m, data = as.data.frame(df_iso_work_dd), plot = F)

pt_DL_mahalanobis_t_4_est <- pt_DL_mahalanobis$estimates[1]
pt_DL_mahalanobis_t_3_est <- pt_DL_mahalanobis$estimates[2]
pt_DL_mahalanobis_t_2_est <- pt_DL_mahalanobis$estimates[3]



pt_DL_mahalanobis_t_4_high <- pt_DL_mahalanobis_t_4_est+(pt_DL_mahalanobis$standard.errors[1]*1.96)
pt_DL_mahalanobis_t_3_high <- pt_DL_mahalanobis_t_3_est+(pt_DL_mahalanobis$standard.errors[2]*1.96)
pt_DL_mahalanobis_t_2_high <- pt_DL_mahalanobis_t_2_est+(pt_DL_mahalanobis$standard.errors[3]*1.96)
pt_DL_mahalanobis_t_4_low <- pt_DL_mahalanobis_t_4_est-(pt_DL_mahalanobis$standard.errors[1]*1.96)
pt_DL_mahalanobis_t_3_low <- pt_DL_mahalanobis_t_3_est-(pt_DL_mahalanobis$standard.errors[2]*1.96)
pt_DL_mahalanobis_t_2_low <- pt_DL_mahalanobis_t_2_est-(pt_DL_mahalanobis$standard.errors[3]*1.96)

Results_DL_mahalanobis <- PanelEstimate(sets = any_closure_binary_control0_m, data = as.data.frame(df_iso_work_dd), se.method = "conditional")
summary(Results_DL_mahalanobis)


# Generate Figure 4

plot_table_NN <- rbind(pt_DL_NN_t_4_est, pt_DL_NN_t_3_est, pt_DL_NN_t_2_est, 0, as.data.frame(Results_DL_NN$estimates))
colnames(plot_table_NN) <- "estimate"
plot_table_NN$high <- c(pt_DL_NN_t_4_high, pt_DL_NN_t_3_high, pt_DL_NN_t_2_high, 0, Results_DL_NN$estimates + (Results_DL_NN$standard.error*1.96))
plot_table_NN$low <- c(pt_DL_NN_t_4_low, pt_DL_NN_t_3_low, pt_DL_NN_t_2_low, 0, Results_DL_NN$estimates - (Results_DL_NN$standard.error*1.96))
plot_table_NN$Matching <- "No Matching" 
rownames(plot_table_NN)[1:5] <- c("t-4","t-3","t-2","t-1","t")
plot_table_NN$Week <- rownames(plot_table_NN)


plot_table_R <- rbind(pt_DL_mahalanobis_t_4_est, pt_DL_mahalanobis_t_3_est, pt_DL_mahalanobis_t_2_est, 0, as.data.frame(Results_DL_mahalanobis$estimates))
colnames(plot_table_R) <- "estimate"
plot_table_R$high <- c(pt_DL_mahalanobis_t_4_high, pt_DL_mahalanobis_t_3_high, pt_DL_mahalanobis_t_2_high, 0, Results_DL_mahalanobis$estimates + (Results_DL_mahalanobis$standard.error*1.96))
plot_table_R$low <- c(pt_DL_mahalanobis_t_4_low, pt_DL_mahalanobis_t_3_low, pt_DL_mahalanobis_t_2_low, 0, Results_DL_mahalanobis$estimates - (Results_DL_mahalanobis$standard.error*1.96))
plot_table_R$Matching <- "Matching + Refinement" 
rownames(plot_table_R)[1:5] <- c("t-4","t-3","t-2","t-1","t")
plot_table_R$Week <- rownames(plot_table_R)

plot_table_DL <- rbind(plot_table_R, plot_table_NN)
plot_table_DL$Week <- factor(plot_table_DL$Week, levels = c("t-4", "t-3", "t-2", "t-1", "t", "t+1", "t+2", "t+3"))

plot_table_DL$`Estimation Strategy` <- factor(plot_table_DL$Matching, levels = c("No Matching", "Matching + Refinement"))


# Plot
pd = position_dodge(.6)   
ggplot(plot_table_DL, aes(x = Week, y = estimate, color = `Estimation Strategy`)) + geom_hline(yintercept = 0, color = "gray40", linetype="dashed") +
  geom_point(shape = 16, size  = 6, position = pd) +
  geom_errorbar(aes(ymin  = low, ymax  = high), width = 0, size  = 1.45, position = pd) +
  theme_bw() + theme() + ylab("Estimated Treatment Effect on Attitudes towards Refugees") + scale_color_manual(values= c("gray70", "#377EB8")) + theme(legend.text=element_text(size=20), axis.text = element_text(size = 20), axis.title = element_text(size = 20), text = element_text(size = 20), plot.title = element_text(hjust = 0.5)) + ggtitle("Border Closures") 

ggsave(path = "~/work/November_2024/R_out/main_plots/", filename = "Refugees_MainPlot.pdf")

